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ABSTRACT 

The large-scale structure in the Universe is believed to arise out of small 
random density perturbations generated in the very early Universe, that 
are amplified by gravity. Large and usually intricate N-body simulations 
are typically employed to model the complex nonlinear dynamics in a self 
gravitating medium. We suggest a very simple model which predicts, on large 
scales, the correct density and velocity distributions. The model does not 
involve an explicit computation of the gravitational force. It is based on a 
simple transformation of the variables and local conservation laws of mass 
and generalized linear momentum. The model demonstrates that the overall 
appearance of large-scale structure in the Universe is explicitly determined by 
the initial velocity field and it reveals the most significant large-scale effects of 
gravity on the formation of structure. 



Subject headings: cosmology: theory - large-scale structure of universe 



1. Introduction 

Understanding the formation of the large-scale structure (LSS) in the Universe is one 
of the most important problems in modern cosmology (see e.g. | ZePdovich fc Novikov 1983 , 
Peebles 1993| ). On scales of galaxies complex nonlinear physical processes like turbulence, 
cooling and heating of the gas, star formation and supernova explosions, determine the 
shape and other properties of objects. However, on the scale of superclusters these 
processes become considerably less important and it is reasonable to assume that primarily 
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gravity shapes the structures. In addition, a variety of observational evidence supports 
the hypothesis that most of the mass in the Universe is in the form of weakly interacting 
particles that feel only gravity. Therefore, a simple physical model assuming that the 
Universe is filled with a cold dust-like gas only interacting gravitationally is most often used 
for studying the LSS formation. The baryonic component plays no crucial role except on 
the scale of galaxies. 

When the amplitude of the density fluctuations is small, their growth is adequately 
described by the linear theory. However, as soon as the amplitude of the perturbations 
reaches the order of unity (a^ = {[{p — ~ 1, where p is the density field and p its 

mean) the linear theory breaks down (see e.g. | Zel'dovich fc INovikov 1983| , [Shandarin fc 



Zel'dovich 1989| ). The most widespread method to deal with the complex dynamics at the 



nonlinear stage is to run N-body simulations generating the initial condition as a realization 
of a Gaussian random process (|Klypin fc Shandarin 1983 ). In N-body simulations of this 



type the gravitational forces generated by the density distribution is calculated at each 
time step. The trajectory of every particle is integrated in a self-consistently varying 
gravitational field. Cosmological N-body simulations have played the most significant role 
in testing (and in most cases rejecting) the models for dark matter. Here we are interested 
in a different aspect of the problem of the LSS formation. We are trying to formulate simple 
physical macroscopic principles controlling the nonlinear stage of gravitational clustering. 

Long ago Zel'dovich (1970) suggested a very elegant, analytical approximation to 
describe the beginning of the nonlinear stage in cosmological scenarios assuming smooth 
initial conditions. Quantitatively it can be expressed as a requirement that the initial 
power spectrum of density fluctuations has a steep cutoff on small scales (steeper than 
P{k) oc A;~^). Mathematically, the Zel'dovich Approximation (ZA) is a one step mapping 
from the Lagrangian space into the Eulerian space at a time t, given by 

r(q,t)=a(t)[q-D(t)V<I'(q)] (1) 

where r and q are the Eulerian and Lagrangian coordinates, respectively; a{t) is the scale 
factor describing the homogeneous expansion of the Universe; D{t) is a known function of 
time describing the growth of perturbations; and $(q) is the potential of the initial velocity 
field: vq oc — V$(q). With the aid of the above mapping one can calculate the density at 
the final time t using the conservation of mass. However, cosmological observations almost 
certainly exclude scenarios having no perturbations on small scales. Small scale power is 
required to explain the existence of quasars and galaxies at very high redshifts. 



Two modifications — the Truncated Zel'dovich Approximation (TZA) ([Kofman et al 



1992| , Poles et al 1993D and the Adhesion Approximation (AA) ( [Gurbatov, Saichev fc 



Shandarin 1989| ) — have been suggested in an attempt to extend the scope of the ZA and to 



- 3 - 



make it more useful for general cosmological scenarios (for a brief review see e.g. Sliandarm 



1994| and for more exhaustive discussion of various approximations see [Satliyaprakash et 



al 1995| , [Saiini fc (Joles 1995| ). A comparison with gravitational N-body simulations shows 



that these two approximations fairly accurately describe nonlinear gravitational clustering 
( [Melott, Shandarin fc Weinberg 1994| , ^athyaprakash et al 1995| ). Here we describe 



another approximation to the nonlinear gravitational evolution of the density and velocity 
perturbations which is numerically almost as simple as the ZA but completely universal and 
much superior to any approximation method suggested so far. It consists of two elements: 

• a transformation of variables and 

• an assumption that in the dense regions the local gravitational interaction can be 
approximated by the diffusion of a generalized momentum (see below for a definition 
of the generalized momentum). 

The latter assumption is similar to the AA but in this model the generalized momentum is 
locally conserved. 



2. conserving Momentum Approximation (COMA) 

It is well known that one can exclude many effects of the homogeneous and isotropic 
expansion of the Universe by introducing the comoving coordinate x = r/a{t) and peculiar 
velocity Vp = u — H(t)r where u is the physical velocity and H{t) = a/a is the Hubble 
parameter characterizing the rate of expansion of the Universe. In addition, it is convenient 
to parameterize time by D{t), which describes the growth of perturbations, and rescale the 
density and peculiar velocity ( purbatov, Saichev fc Shandarin 1989|) by 

p(x, t) = a-^r]{^, D{t)), vp(x, t) = (al))v(x, D{t)). (2) 

In terms of the new variables the dynamics of gravitational clustering is described by (i) 
the conservation of mass, (ii) a dynamical equation analogous to the Euler equation 

dvj dvj _ 3 VLq ( ^V , \ (o\ 

where Qq = SuGp/SHq is the dimensionless mean density of the Universe, and (iii) an 
equation for the gravitational field ip analogous to the Poisson equation (for derivations see 
e.g. Shandarin 1994| ). 
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As in the ZA the COMA assumes that the gravitational force is approximated by the 
velocity 

dif/dxiK. -Vi{^,D) (4) 

which sets the right hand side of eq.|0| to zero. Thus, on large scales the model mimics the 
dynamics of a self-gravitating medium by the effectively kinematic model of non-interacting 
medium. 

However, on small scales the approximation assumes that particles inelastically collide 
resulting in the exchange and diffusion of the generalized momentum defined as 

P = '7V = -^pvp, (5) 



which mimics the effects of gravity in dense regions ([Shandarin fc Zel'dovich 1989|). 



In the AA the right hand side of eq.||^ is supplemented by the viscosity term resulting 
in Burgers' equation 

dVi dvi 2 fa\ 

In the AA u is assumed to be constant which results in the violation of the momentum 
conservation ( [Kofman fc Shandarin 1990| ) and in this respect the model is not physical. The 
COMA assumes physical viscosity described by the equation similar to the Navier-Stokes 
equation (here we ignore the 2nd viscosity) 

dD ^ dxk 7] dxk dxk dxi 3 dxi 

where /x is the dynamical viscosity. The Navier-Stokes equation conserves the linear 
momentum automatically but does not have an analytical solution in the general case; also 
it may generate vorticity. 

According to the ZA in the single stream flow regime each particle moves along straight 
line with constant velocity. However, as soon as a particle enters a multi-stream flow 
region its trajectory becomes very complex, resembling a random walk. Collectively this 
type of motion can be labeled as a sort of violent relaxation ( [Lynden-Bell 19(171 ) and can 



be roughly approximated by the diffusion of the momentum. It is convenient to split the 
force exerting on the particle into two parts: The first component represents its interaction 
with the other particles of the clump of which it is a member and the second represents its 
interaction with particles that belong to other clumps or voids. One can reasonably assume 
that the irregular component of the motion is mostly determined by the gravitational force 
induced by the particles from the same clump. If so, the momentum of the clump must be 
approximately conserved in the interaction of this type. The change of the momentum of 
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the clump is due to interaction with the surrounding matter that can be roughly described 
by eq. [4]. In practice, the diffusion of momentum is important in all regions of high density 
(filaments and pancakes) and aids in preventing the formation of multi-stream flows as 
in the AA. The major physical difference between the COMA and AA consists in that 
the COMA conserves momentum locally and the AA, based on Burgers' equations, does 
not (instead, the latter conserves the velocity). We believe that this an advantage of the 
COMA. 

Another advantage of the COMA is mainly practical. The numerical code realizing 
the model is extremely simple and efficient. It operates with the density and velocity 
distributions (no particles!) on a cubic grid in the Eulerian space: r] = ri{'x,D) and 
V = v(x, D). At each time step it calculates the flows of mass and momentum from a mesh 
cell to its neighbors and computes the new density and velocity fields on the grid using a 
cloud-in-cell algorithm (see e.g. Hockney fc Eastwood 1988 ) thereby explicitly conserving 



the mass and the generalized momentum (simulation of perfectly inelastic collisions). Thus, 
one circumvents the problem of having to compute the gravitational force after each time 
step. A full description of the code as well as the code itself will be published elsewhere; 
here we list its major features. The algorithm implementing the new approximation is: 



• universal in terms of the initial and boundary conditions as well as the shape of the 
computational box, provided that the initial velocity field is somehow generated; 

• extremely simple both conceptually and practically; 

• very economical in terms of memory (just four arrays, each as big as the size of the 
box, are sufficient to implement the algorithm); 

• very fast; and 

• local, therefore 100% parallelizahle. 



Generalization of the algorithm to deal with particles, rather than densities on the grid, 
as well as introducing partially inelastic collisions is straightforward but would require 
additional memory and slightly slow down the code. 

The only disadvantage of the COMA (compared to the AA) is that it does not have an 
analytic solution. 
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3. Comparison with N-body simulations 

In order to test the strength of our approximate model we compare it with gravitational 
N-body simulations. The N-body simulations are performed using 128'^ particles on a 
128^ mesh with periodic boundary conditions. The power spectrum of primordial density 
fluctuations is assumed to be a simple power law (|5fcp oc fc", where 5k is the Fourier 
transform of the density contrast 5) covering the range n = -|-1, 0, — 1 and —2. The initial 
fluctuations are evolved gravitationally and comparisons are made when different scales, 
determined by the parameter k^i (defined by the equation cr| = a? Jq"' P{k) d?k = 1), go 
nonlinear: k^i = 64, 32, 16, 8 and 4 (in units of the fundamental frequency kf = 271 /L, where 
L is the size of the box). For a detailed discussion of the N-body simulations see [Melott & 



[Shandarin 1993| . At each epoch the density field is computed on a reduced grid of size 64^ 



by a cloud-in-cell method. 

We run the model code with identical initial conditions and compare its results with 
those of N-body simulations at different stages. Here we illustrate this comparison for three 
different values of the power-law index of the spectrum of fluctuations, namely, n = 0, —1 
and —2, and at two stages kni = 8 and 4. Although it is unlikely that such power-law models 
can explain the real universe, they serve as toy models to understand generic features of 
gravitational instability. The following normalizations provide a rough idea of the scales 
involved and how such toy models may relate to the real Universe. One can view the 
density fields corresponding to different stages in the evolution of such power law density 
fluctuations as equivalent to the density distribution at the present epoch but obtained after 
smoothing with a top-hat filter of progressively smaller radii Rth, with smaller smoothing 
lengths corresponding to later epochs. For instance, the epochs kni = 8 and 4 correspond to 
smoothing radii Rth ~ 2, and 1 Mpc, respectively, within volumes (200 Mpc)^, 
and (100 /t-i Mpcf. 

In Fig. 1 we compare the N-body density distributions (open circles) with those 
obtained using the COMA (crosses) in a typical two-dimensional slice, one mesh size thick, 
in the n = —1 case, when the scale of nonlinearity corresponds to kni = 8. The panels show 
sites of density contrast (6) larger than a certain threshold {6c). A filled square is a site 
where both the N-body and model densities are larger than the threshold. One can see 
that there is generally a good overall agreement between the two distributions at nonlinear 
but not extremely high density thresholds. However, the COMA shows a rather smoother 
distribution and the density peaks that it produces are somewhat lower, especially at 
higher densities. This behavior is not unexpected: local nonlinear gravitational effects, 
ignored by the COMA, make the density distribution more clumpy and the clumps are 
more compact. A detailed comparison shows that the vast majority of density peaks in 
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Fig. 1. — Sites with densities exceeding p > 1,2,4, and 8 times the mean density in the 
n — —1 model at the kni — 8 stage, are shown. Empty circles correspond to sites where only 
the N-body densities are larger than the density threshold, crosses represent sites where only 
the COMA densities are larger than the threshold and filled circles show sites where both 
the distributions are larger than the threshold. 

the N-body simulations have their counterparts in the COMA often shifted by a couple 
of mesh sites and of slightly lower amplitudes. In order to quantify the similarities and 
discrepancies of the two density distributions in Fig. 2 we plot the power spectra (left 
panels) and probability distribution functions (right panels) of the N-body (solid lines) 
and the COMA (dashed line) density fields, for three values of the power-law index n — 0, 
— 1 and —2, (panels from top to bottom) and two epochs kni — 8 (bottom pair of curves) 
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Fig. 2. — The left column shows the power spectrum in three models at two stages each: 
the bottom curves correspond to the kni — 8 and the top curves correspond to the kni — 4. 
Solid lines correspond to the N-body simulations and dashed hues to the COMA. Similarly, 
the right column shows the density distribution functions normalized to the total number of 
mesh sites (64^). In order to avoid overlapping the curves for the kni — 4 stage are multiplied 
by a factor of 10. 

and 4 (top pair of curves). Comparisons of the power spectra show that the COMA does 
not have enough power on small scales which is obviously related to a relative smoothness 
of the density fields in the COMA. Similarly, the density distribution functions show that 
when kni — 8; sites with densities above 10 (in units of the mean density) are less abundant 
in the COMA. However, at a later stage when kni — 4 the COMA density distribution 
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functions reproduce the N-body distribution correctly up to densities p = 30, but at larger 
densities (not shown) they again fall down faster than the ones in the N-body simulations. 
A similar comparison of N-body with other successful models, such as the AA or the TZA, 
does not reveal as great an agreement as between COMA and N-body ( Melott, Shandarm 
fc WeinbergT994D. 



4. Summary 

Summarizing, we conclude that the proposed approximation mimics the large-scale 
gravitational evolution at late nonlinear stages quite well (positions of clumps, filaments, 
pancakes) - better than any other known approximation. This implies that the two simple 
assumptions: 

• the generalized gravitational force on large scales (/ > k~l) equals the velocity : 
dip/dxi ~ — t>j(x, D) and 

• the conservation of mass and the local {I < k~l) conservation of the generalized 
momentum, p = r^v. 



explain fairly well the nonlinear gravitational clustering on large scales. 

In a very general sense the COMA falls in the class of sticky particle methods 
used in the numerical hydrodynamics and also resembles the lattice gas models used in 
modeling turbulence and similar phenomena. An interesting question is whether the COMA 
guarantees complete hierarchical clustering or not. Due to numerical errors some small 
clumps may miss merging with the larger ones passing by without the collision. We have 
not seen this phenomena in our simulations and believe that it should not be a serious 
problem for the method. 

We believe that after a thorough testing the approximation can be a practical tool for 
cosmological studies of large-scale processes which do not require a resolution better than 
a few Mpc, such as large-scale streaming velocities, spatial distribution of rich clusters of 
galaxies, statistics of voids, etc. In such low resolution calculations the COMA can be more 
efficient than N-body simulations if very large volumes and large statistical ensembles are 
required. 
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